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Abstract-In this article a new approach is considered for 
implementing multi-grid methods in iterative operator splitting 
methods for differential equations. The underlying idea is to 
embed fast multi-grid methods to accelerate the iterative splitting 
schemes. We concentrate on convection-diffusion-reaction 
equations and treat a multi-scale problem. The equations are 
spatially discretised with finite difference or finite volume 
methods. The main problem of iterative solvers is their neglect of 
multi-scale physics: different scales are not resolved on different 
spatial scales. Here we apply multi-grid methods to obtain the 
resolution of the optimal scale, which can be solved by an 
iterative splitting scheme. We discuss the embedding of such 
spatial- and time-scale methods and their underlying analysis. 
Numerical examples are given to verify the numerical schemes. 
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I. INTRODUCTION 

Our motivation is to solve multi-scale problems that occur, 
for example, in transport-reaction problems in porous media. 

In the last years, the interest in numerical simulations of 
multi-scale problems that can be used to model different 
physical behaviors, e.g., potential damage events, has 
significantly increased, but the adaptation of this methodology 
to such problems is delicate, see [39, 40]. 

We will concentrate on simplified models and give an 
account of the ideas needed to adapt the underlying splitting 
schemes, see [15]. While these splitting schemes worked well 
for single scale problems, we have to embed multi-grid ideas 
or multi-step ideas into these schemes to extend them to be 
multi-scale solvers. We will deal with iterative splitting 
methods and how to decompose such scale problems and 
accelerate the convergence rates with multi-scale methods, so 
that we can overcome such scaling problems. 

Moreover the embedding part is important: while 
dealing with multi-grid methods, we have to extend 
the underlying splitting analysis. We will discuss the 
needed algorithmic ideas. 

The novelty is to cheaply embed multi-grid and multi-step 
methods and apply only cheap iterative schemes to achieve 
higher order results. 

In general, we discuss an elegant way of embedding and 
survey recent methods for multi-scale problems with respect 
to iterative splitting schemes. 

In the following, we describe our model problem. The 



model equation for the multi-scale equations are given as 
coupled partial differential equations. 

We concentrate on a far-field model for a plasma reactor, 
see [16] and [30], and assume a continuum flow and that the 
transport equations can be treated with a 
convection-diffusion-reaction equation, due to a constant 
velocity field: 



- + V-Fc= 0,innx [0,t] 
at L ' J 

F = v - DV, 

c(x,t) = c (x),onfl, 

3c(X,t) 



flu 



= 0, on da x [0, t] 



(1) 

(2) 
(3) 



where c is the particle density of the ionised species, F is 
the flux of the species, v is the flux velocity through the 
chamber and porous substrate, which is influenced by the 
electric field, and D is the diffusion matrix. The initial value is 
given as c and we assume a Neumann boundary condition. 

The aim of this paper is to present a novel iterative 
splitting method that embeds multi-scale methods for spatially 
discretised transport-reaction equations. 

This paper is organized as follows: In Section II, we 
present the underlying splitting methods. The stability analysis 
of the embedded multi-grid schemes with the iterative 
schemes is given in Section III The assembling and the 
algorithms for the embedded multi-grid method are given in 
Section IV. Numerical verifications are given in Section V. In 
Section VI we briefly summarize our results. 

II. ITERATIVE SPLITTING METHODS 

Iterative splitting methods were developed in the early 90s, 
see [28, 41]. 

The idea is to decouple the equations into two or more 
equations to save computational time, without changing the 
previous results. We obtain inhomogeneous partial differential 
equations and solve them with the appropriate methods. 

We consider the following differential equation problem, 
where the operators A and B are spatially discretised: 



da 

at 



= Au + Bu 



(4) 



where the initial conditions are un = u(tn). The 
solutions correspond via the spatial discretisation by u(t) 

= (ul (t), . . . , urn (t))T = (c(xl , t), . . . , c(xm , t)T and m 
is the number of spatial discretisation points. The operators 



JMMF Vol.1 Issu.2 2012 PP.1-10 www.sjmmf.org ©2012 World Academic Publishing 

1 



Journal of Modern Mathematics Frontier 



JMMF 



A and B are matrices corresponding to the discretised 
convection and diffusion operators in (1) and their boundary 
conditions. Hence, they can be considered as bounded 
operators with a sufficient large spatial step Ax > 0. 

In the following we discuss the different methods based on 
the following definition: 

Definition 1. We define a sequential method as a serial 
method, while each step can only be done successively (e.g., 
each processor has to wait till the results of the previous 
computation are done). We define a parallel method as a 
multitasking method, which can be executed one piece at a 
time on many different processing devices, and then put back 
together again at the end to get the correct result, (e.g., each 
processor can independently compute a result). 

A. The Sequential Iterative Splitting Method 

The classical sequential operator splitting methods, known 
as Lie-splitting or Strang-Marchuk splitting, have several 
drawbacks besides their benefits, see [38, 31]. For instance, 
for non -commuting operators there might be a very large 
constant in the splitting error which requires the use of an 
unrealistically small time step. Also, splitting the original 
problem into different sub-problems with one operator, i.e., 
neglecting the other components, is physically questionable. 

In order to avoid these problems, one can use the iterative 
operator splitting method on the interval [0, T]. This algorithm 
is based on an iteration with a fixed splitting discretisation 

step-size x . On every time interval [t n , t n ] the method 
solves the following sub-problems consecutively for i = 0, 

2, . . . 2m. 



auj(t) 
at 

du i+1 (t) 

at 



= AUi(t) + BUi.iOO.withUiO: 11 ) = u n ,t 6 (t n ,t n+1 ), (5) 



= AUj(t) + Bu i+1 (t), with u i+1 (t n ) = u n , t 6 
(t n ,t n+1 ), 



(6) 



where u n is the known split approximation at the time level t 

= t n (see [9]).. This algorithm is an iterative method which 
involves in each step both operators A and B. Hence, there is 
no real separation of the different physical processes in these 
equations. 

Remark 1. For the presented iterative splitting method, we 
have a serial algorithm and we cannot use parallel methods in 
this version. For the sake of efficiency, we will modify this 
method to enable parallelisation. 

B. The Parallel Iterative Splitting Method 

To parallelise the standard iterative splitting method, see 
Section II- A, we propose a decoupled version. 

The iterative scheme can be done in a Jacobian form, so 
that the two parts of the algorithm can be computed 
independently. 

In the next iteration steps, we use the previous solutions 
and improve them. 

We apply the parallel iterative splitting method for i = 0, 

2, . . . 2m, and we have: 

^- = Au;(t) + BUi.iOO.with Uj(t n ) = u n , t 6 (t n ,t n+1 ) (7) 
u (t n ) = u n ,u i _ 1 (t) = 0, 



flUi + l(t) 

at 



= A Ui _ 2 (t) + Bu i+1 (t),withu i+1 (t n ) = u n ,t 6 



(8) 



(t n ,t n+1 ), 

u_ 2 (t) = 

where u n is the known split approximation at the time level 
t=tn. This algorithm is an iterative method where each step 
involves both Operators A and B. Hence, there is no real 
separation of the different physical processes. 
The first iterative steps are: 

3u (t) 



at 



■ = Au (t),withu (t n ) = u , 



^^ = B Ul (t), 

at 1V J 

withu^t") = u", 



(9) 
(10) 
(11) 



and 



au 2 (t) 
at 

3u 3 (t) 



= Au 2 (t) + Bu 1 (t),withu 2 (t n ) = u , 
= Au (t) + Bu 3 (t), 



at 

withu 3 (t n ) = u n 



(12) 



Remark 2. The effect of the parallelisation is to obtain the 
accuracy of the sequential iterative splitting method with one 
more iterative step. 

C. The Multi-grid Algorithm 

To describe the multi-grid method we first implement the 
two-grid method. Afterwards, we extend it recursively to the 
multi-grid method. The smoother grid Level 1 is denoted by 
SI. The two-grid method is defined as follows: 



Mf G := S t V2 Mf GG S^ 1 



(13) 



where vl denotes the pre-smoothing steps and v2 denotes 
the post-smoothing steps. The correctionMf 00 on the coarse 
grid is defined by: 



Mf GG : = 



1 - pA^rA 



(14) 



The passage to the multi-grid method is done in such a 
way that Al-1 of the coarse grid in Equation (14) is not 
inverted exactly, but the two-grid method is invoked y-times 
to solve the equation systems at the grid level 1-1. The 
system of equations is only solved on the coarsest grid. 



The multi-grid method is defined as: 
Mg 10 := 0, 
M, MG := Mf G 
m mg : =Ml ZG +S 1 V2 p(M 1 M . G ) Y Ar_ 1 1 rA 1 S 1 Vl 



(15) 
(16) 
(17) 



where vl denotes the pre-smoothing steps and v2 the 
post-smoothing steps. The corrections Mf cc of the coarse 
grid are: 

m mgg :=1 _ p Q _ ^M^_ G l y)Al} 1 YA l . (18) 

This approach is called the multi-grid method. For a 
choice of y =1 one speaks of a V-cycle, for y = 2 it is a 
W-cycle. 
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The multi-grid algorithm is given by 
Algorithm 1 Linear multi-grid cycle 

MG (x h b lt l) 

{ 

if {I == 0) 

{ 

x = A^bg, exact solving on coarse grid. 

} 
else 

{ 

x x = S Vl {x b bi); v x -pre — smoothing steps. 
b t _ 1 = rb t ; defect restricted on 

next coarser grid. 
x h = 0; 
for (i=0; i < y; i++) 

{ 

Ci-i = 0; 

MG(ci_-y, bi-i, I — 1); y — fold recursive 
invoking of correction of coarse grid 

x l-l = x l-l + C|-1j 

Vi = bi_i-A,_ 1 c 1 _ 1 . 

} 

q = pxi_!; interpolaton of correction 

of the next coarse grid, 
x, = x, + c,; 
b, = b,-A,c,; 
X[ = S V2 (X[,b[); v 2 post — smoothing steps. 

} 

} 

The refinement is done using the strategy of [18]. By 
assuming a linear effort for the smoothers as well as for the 
grid-transfer operators, one obtains a linear effort for 
Algorithms 1, if y < 3, cf. [18]. The proofs of convergence for 
the W-cycle were given in [17]. The topic of convergence will 
not be further discussed, for an overview we refer to [44]. 

Subsequently, the multi-grid cycles are illustrated in Fig. 1. 
For a further treatment of the topic of multi-grid methods, we 
refer to [2, 17, 18]. 




- 

A 



-L Level ■ 

V-Zyklus 



A/W 



• * Voiglatten 
• Nacliglatten 
■ ExaktLosen 



W-Zyklus 



Fig. 1 Multi-grid cycles: V- and W-cycle 

D. Iterative Splitting Method with Embedded Multi-grid 
Method 

The following algorithm is based on embedding the 
multi-grid method into the operator splitting method. The 
iteration with fixed splitting discretisation step- size x . On the 
time interval [tn , tn+1 ] we solve the following sub-problems 
consecutively for i = 0, 2, . . . 2m. (cf. [24] and [9].) 



^P = ACi(t) + P 1a " 1b + Bq_ 1 (t),withq(t n ) = c n , (19) 

££i±iW = r1a-1b Ac^t) + Bc i+1 (t), with c i+1 (t n ) = c n , (20) 

where c (t n ) = c n , c—] = and c n is the known split 

approximation at t = t n . We assume A is the fine spatial 
discretised operator on level IA , where B is the coarse 
discretised operator on level IB . 

The operators are coupled by the restriction and 
prolongation operators: 



A = P'a-'ba 

■"coarse L **' 

B fme =P 1 A-iB B , 



(21) 
(22) 



where R is the restriction and P the prolongation operator. 

Algorithm 2 We assume two spatial operators A, B, which 
are discretised by finite difference of finite element methods. 

We solve the time-discretisation of our equations 

^ = Aq(t) + P lA - 1 BBq_ 1 (t),withq(t n ) = c n , (23) 

££i±l^ = rIa-ib Aci(t) + Bc i+1 (t), with c i+1 (t n ) = c n , (24) 

where CO (t n ) = c n , c—] = and c n and cn is the known 
split approximation at time-level t = t n by integration with 
respect to time and obtain 
(I - A)Ci(t n+1 ) = Ci(t n ) + P'A-'BBCi.iCt 1 ^ 1 ), 

withq(t n ) = c n (25) 

(I - B)c i+1 (t) = c i+1 (t") + R 1 A-'BAc 1 (t n+1 ), 

withq +1 (t n ) = c n , (26) 

We have the multi-grid equations: 
L iCi (t n+1 ) = q(t) n + P 1 A- 1 BBq_ 1 (t n+1 ), 

withq(t n ) = c n , (27) 

L i+1 q +1 (t) = q +1 (t") + R'A-iBAqCt"- 1 " 1 ), 

withq +1 (t n ) =c n , (28) 

E. The Waveform Relaxation Method 

A further method to solve large coupled differential 
equations is the waveform relaxation scheme. 

The iterative method was discussed in [41] and can be 
done either with Gaussian or Jacobian form. 

We deal with the following ordinary differential equation 
or assume a semi -discretised partial differential equation: 

ut =f(u, t), in (0,T), u(0) = vO 
where u = (ul , . . . , um)t and / (u, t) = (fl (u, t), . . . , fm (u, 

t)f. 

We apply the waveform relaxation method for i = 0, 1, . . . m 

and have: 

du 1;i (x,t) 



at 



du 2 ,i(x,t) 
dt 



— fl( u l,i< u 2,i-l< ■■■'U mi _ 1 J 

withu^Ct") = %(t n ) 

= f2(. u l,i-l< u 2,i< u 3,i-l< ■■■' U m,i-lJ 

withu 2J (t n ) = u 2 (t n ) 



(29) 



(30) 
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du mi (x,t) 

at 



= f2(ui,i-l, 



-l,i-l- u m,iJ 



withu m4 (t n ) = u m ,(t n ) 



(31) 



where for the initialisation of the first step we put u 1 _ 7 (t) 

= u,(t n ),..., u m ,_ 1 (t) = ujt n ). 

We reduce to two equations and reformulate the method to 
our iterative splitting methods. 



So we consider 



da 



^ = f 11 (u 1 ,t)+f 12 (u 2 ,t),in(0,T) 



= f 21 (u 1 ,t)+f 22 (u 2 ,t),in(0,T) 



5u2 

at 



u(0) = v 



where u = (u\ , u2 ) ■ 
Our notation for the operator equation is 
^ = A(u) + B(u),in(0,T) 
u(0) = v 

Where 

■*■> = te) 

The iterative splitting method as waveform relaxation 
method written for 

i = 0, 1, . . . m as: 

7f = A(u li ,u 2i _ 1 ) + B(u li _ 1 ,u 2i ) 

with u 14 (t n ) = Ul (t n ) andu 24 (t n ) = u 2 (t n ) 



(32) 

(33) 
(34) 



(35) 
(36) 

(37) 
(38) 



(39) 



where for the initialisation of the first step we put u } _ 1 (t) = 

u,(t n ),u 2 _ J (t) = u 2 (t n ). 

Convergence analysis In the following we formulate the 
iterative splitting method as a waveform relaxation method 
and apply the proof techniques of those relaxation schemes. 

The waveform relaxation method is given by 
diij 



dt 



= PU; + Qu^ + f, 



u,(t n ) = u(t n ), 



(40) 
(41) 



u i+1 (t n ) = u(t"), 



(45) 



where P, Q are matrices given by spatial discretisation, e.g., 
P is the convection part of Q the diffusion part. 

But we can also make an abstract decomposition, taking 
into account that A = P + Q, where P is the matrix with small 
eigenvalues and Q is the matrix with large eigenvalues. 

Theorem 1 We have the iterative operator splitting 
methods given as the following outer and inner iterations. 



Outer Iteration: 



[ - = PU; + QU;.! + F, 



dU 
dt 



u,(t n ) = u(t n ), 



(46) 
(47) 



where P and Q are the diagonal and outer-diagonal 
matrices of the splitting methods. 



Inner Iteration: 



duj 

dt 



= PUj + Qu H1 + f, 



Uj (t n )=u(t n ), j = l J 

^ = Pu, + Qu k + f, 
u k (t n ) = u(t n ), k=l 



(48) 
(49) 
(50) 
(51) 



We have a convergent scheme for K bounded in a 
Banach-space: 

p(K) := lim ||K K || 1/K (52) 



where p(K) < 1 is the spectral radius of K and is given by 
the variation of constants of P and Q. 



Proof. The outer iteration is given by 

dUj 



dt 



= PU; + QU;.! + F, 



UiCt") = U(t n ), 



(53) 
(54) 



where P and Q are the diagonal and outer-diagonal 
matrices of the splitting methods, which can be solved by the 
variation of constants. 



We introduce the linear integral operator 
KU(t) = J t exp((t - s)P) QU(s)ds, 



(55) 



where A = P + Q, e.g., P is the diagonal part of A (Jacobi 
method). 

Here the splitting method is done abstractly with respect to 
the Matrix A. The method is considered an effective solver 
method with respect to the underlying matrices. 

The reformulation of the iterative operator splitting 
method is given by 



^ = Pu ; + QUi.! + f, 

Ui (t n ) = u(t n ), 

du i+1 



dt 



= PU; + QU i + 1 + f, 



(42) 
(43) 
(44) 



And we have 

Uj = KU;.! + 0, 



(56) 



= exp(tP)U(t n ) + J t exp((t - s)P) F(s)ds, (57) 

For the convergence it is sufficient to show that K is 
bounded in a Banach-space: 

p(K) := lim ||K K || 1/K (58) 

K — >co 

where p(K) < 1 is the spectral radius of K 

This is given by the bounded operators, see [33, 34]. 
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III. ERROR ANALYSIS: STABILITY THEORY FOR THE ITERATIVE 
SPLITTING METHOD WITH MULTI-GRID METHODS 

The following algorithm is based on embedding the 
multi-grid method in the operator splitting method. The 
iteration with fixed splitting discretisation is step-size x. On 

the time interval [t n , t n+ ] we solve the following 
sub-problems consecutively for i = 0, 2, . . . 2m. (cf. [24, 9].) 

^-^ = Aq(t) + P'A-'BBq.iCO.withqCt 11 ) = c n , (59) 

= R'a-'b ACi(t) + Bc i+1 (t), with c i+1 (t n ) = c n , (60) 



at 
3q+i(t) 



at 



where c (t n ) = c n , c_ } = and c n is the known split 

approximation at time-level t = t n . We assume A to be the 
fine spatial discretised operator on level IA and B to be the 
coarse discretised operator on level IB. 

The operators are coupled by the restriction and 
prolongation operators: 



A = R 1a_1b A 

■"coarse lv -"' 

Bflne =P lA " lB B, 



(61) 
(62) 



where R is the restriction and P the prolongation operator. 

Theorem 2 Let us consider the abstract Cauchy problem 
in a Banach space X 

d t c(t) = Ac(t) + P'a-'bBcOO, < t < T 

c(0) = c (63) 

where A, P " ~' B B, A + P lA ~' B B : X -> X are given 
linear operators being generators of the Co-semigroup and c 
E Xisa given element. Then the iteration process (59)-(60) 
is convergent and the rate of convergence is of higher order. 

Proof. We assume A + P lA ~ lB B E L(X) and assume a 
generator of a uniformly continuous semi -group, hence the 
problem (63) has a unique solution c(t)=exp((A + P A ~ lB 
B)t)c . 

Let us consider the Iteration (59)-(60) on the sub-interval 
[tn, tn+1]. For the local error function ei(t) = c(t) - ci (t) we 
have the relations 

a t ej(t) = Aej(t) + P'A-'BBej.iCt), t 6 (t n ,t n+1 ] 

ei (t n )=0, (64) 

d t e i+ i(t) = R'A-'BAqCt) + Be i+1 (t), t 6 (t n ,t n+1 ] 

e i+1 (t n ) = 0, (65) 



And 



for m =0, 2, 4, . . ., with e0(0) = and e-l(t) = c(t). In the 
following we use the notations X2 for the product space X x 
X endowed with the norm 

||(u,v)|| = max{||u|U|v||}(u,v 6 X) The elements 

£j (t)< F; (t) 6 X 2 and the linear operator 

A : X 2 — > X 2 are defined as follows: 



Ei(t) = 



e ; (t) 
-i+i(t) 



,F,(t) = 



P'A-'BBej.iCt) 




A - Ir'a-'ba bJ 



(66) 



Then using the Notation (66), the Relations (64) and (65) can 
be written in the form 

3 t £i(t) = A£i(t) + Fi(t),t 6 (t n ,t n+1 ], Ei (t n ) = 0. 

(67) 
Due to our assumptions, A is a generator of the one-parameter 
Co-semi-group (exp At) t > 0, hence using the variations of 
constants formula, the solution to the abstract Cauchy problem 
(67) with homogeneous initial condition can be written as 

Ei(t) = J exp(A(t- s))Fi(s)ds,t 6 [t n ,t n+1 ] 



(See, e.g., [8]) Hence, using the notation 

HeJL = SUPtg^ntn+ljIlEiOOH 

We have 

l|£ i ll(t)<l|F i |L/ t t n ||exp(A(t-s))||ds = 
llB||||e i _ 1 ||/ t t n ||exp(A(t-s))||ds,te[t",t" +1 ]. 



(68) 
(69) 

(70) 



Since (A(t))t > is a semi-group, the so called growth 
estimation 

||exp(At)|| < Kexp(wt),t > 0, (71) 



holds with some numbers K > and co £ IR, cf. [8]. 

- Assume that (A(t))t > is a bounded or an 
exponentially stable semi-group, i.e., (71) holds with some co 
< 0. Then obviously the estimate 



||exp(At)|| < K,t> 0, 



(72) 



holds, and hence, on the basis of (70), we have the relation 

IM(t) < KlIP'A-'BBllTjIe^lU 6 [t",t n+1 ] (73) 

- Assume that (exp At)t > has exponential growth 
with some co > 0. Using (71), we have 

J t t n ||exp(A(t- s))||ds < K w (t),t 6 [t",t n+1 ] (74) 



Where 



K w (t) = - (exp(w(t - t 11 )) - l), t 6 [t n , t n+1 ] (75) 



Hence 



K w (t) < - (exp(wT n ) - 1) = K Tn + 0(T n 2 ). (76) 

W 11 

The Estimations (73) and (76) result in 

llqlloc =K||P 1 a-'b ||||B||T n ||e 1 _ 1 || + 0(T n 2 ) (77) 

Taking into account the definition of Ei and of the norm k • koo, 
we obtain 



|e,|| = KlIP'A-'BllHBllTjIei.J + 0(T n 2 ) (78) 



And hence 
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||e i+1 || = KiT^He^JI + 0(T n 3 ), (79) 

which proves our statement. 

Operator-splitting method with embedded Jacobian 
Newton iterative method: Newton's method is used to solve 
the nonlinear parts of the iterative operator-splitting method, 
see the linearisation techniques in [24, 25]. We apply the 
iterative operator-splitting method and obtain 

FiCuj) = d t Uj - A(Ui)Ui - BCuj.^Ui.! = 0, 

withu;(t n ) = c n , 

F 2 (u i+ i) = d t u i+1 - ACuJUi - B(u i+1 )u i+1 = 0, 
withu i+1 (t n ) = c n , 

where the time step isr = t n+ — t n . The iterations 
are i = 1, 3, . . . , 2m + 1. c (t) = is the starting solution and 
cn is the known split approximation at time level t = f. The 

results of the methods are c(t n + * ) = u2m + 2 (t n+ .The 
splitting method with embedded Newton's method is given by 

U ( k+ D = u 00 _ D (f^OO))- 1 ( dtU <M - A(u«)u« 

-B(u i «)u«),withD(F 1 (u«)) 

= -(a(u«) + 



aA(u, (k) ) 



,(k) 



a u 



(k) "i 



andk = 0,1,2, ...,K, with Ui(t n ) = c 11 , 
u^ = u« 1 -D(F 2 (u« 1 ))- 1 ( at u« 1 -A(u«)u« 



B(u«)u«)c«), 



withD(F 2 (u^)) = -(B(u^) + ^^u« 1 | 

andl = 0,l,2,...,L,withUi + l(t n ) = c 11 . 

Remark 3. For the iterative operator-splitting method with 
Newton's method we have two iteration procedures. The first 
iteration is Newton's method for computing the solution of the 
nonlinear equations; the second iteration is the iterative 
splitting method, which computes the resulting solution of the 
coupled equation systems. The embedded method is used for 
strong nonlinearities. 

IV. ASSEMBLING OF THE SPLITTING METHODS WITH 
THE EMBEDDED MULTI-GRID METHOD 

In the following we discuss the assembling of our methods. 

A. Crank-Nicolson and Two-grid Method 

In the following we apply the iterative splitting method, 
with a Crank-Nicolson method in time, a second order finite 
difference method in space, and an underlying two-grid 
method as solver. 

We apply the Crank-Nicolson scheme separately to the 
two equations obtained by the two steps of the iterative 
operator splitting method after the discretisation of the 2D 
heat equation. For the ith step of the iterative operator splitting 
method, as interpolation operator we choose the bilinear 
interpolation given by 

v5- 9 - = v 2h 

v 2i,2j v ij 



V2 h i + l,2 j =^K+V i 2 + h 1 , j ). 



-,2j - 2 V2j 
1 



v5-,-,., = -fv 2h + v 2i \,) 

v 2i,2j + l o ^ v ij T y i,]+U 



v 2 h i +1 ,2 j+1 = 4 (vf + vfAj + < +1 + v? + Vi). 

0<i,j.<^-l. 

where the superscript h is associated with the fine grid and 2h 
with the coarse grid. 



Step i: 

,,1+1 _ ,,n i 

u i.j.k u ij,k _ - 1 - 



At 



,n+l 
-M-Uk 



9,,n+l + ,,n+l 
^ u ij,k ^ u i+lj,k 



+ PD, 



,n+l 
"M-lj-l.k-l 



h 2 

"x 

- 2u n+1 



+ u 



i— l,j,k— 1 ^ u i-lj+l,k-l 



h 2 



+ Di 



u i-i,i,k ^Uj j k + u i+1 j 



+ PD 



h 2 

"x 
u i-l,j-l,k-l — 2Ui_i j k _ 1 + U;, 



+l,k-l 



h 2 

Hy 



_£ix2,,n+l , E2_pfi2 n+1 , £ix2 n , £2_ps:2 n 

~ 2 X u ij,k 2 y i-lj.k-1 + 2 x iJ' k 2 Y i-lj,k-l 



Pu 



-D 2 At 



DiAt\ n+1 



+ u 



n+1 



,n+i 1 



_i_ ,,n+i — |_ Pnii+i 

^ ij.k u2 ^ ru i-lj,k-l 



D 2 At 
h 2 



~ ru i-l,j+l,k-l 



-D 2 At 
2h y , 



n + 1 (~ D l At \ 

+ ur+1 J' k V^hTJ 



= 2 ( D i 6 x u U,k + D 2 PS y un_ Uk _ 1 ) + <j, k 



u 



+ u 



n + 1 

i— l,j— l,k 

DiAt 1 



-D,At 



2h 2 



+ 



u r-uk( : 



-DiAt 
2h 2 



+ u 



n + 1 

ij,k 



n+1 



•J* h Y 2 2 



+ 9 l u i-lj-l,k + u i-l,j+l,kJ 



D 2 At 
h 2 

Hy 



1 

4 l 



J.^f,,n+1 4-i,ii+l , n + 1 . n+1 \l 2_ 

r A l u i-lj+l,k T u i,j+l,k T u i-lj+2,k ^ u i,j+2,kJ \ 2h 2 



D,At 



,n+l 



DiAt\ 



4- ll n_t 

+ Ui+1 J' k V 2h 2 / 

-"""n^ 2 ,," 4- 1 nP ^-Li-i-k-i ~ 2u "-ij,k-i + u "-ij+i,k-i 

- 2 U l°x u i,j,k + nU 2 f- 



\i,k - 2 ' 



+ u "i,k 



h 2 

Hy 



D,At\ ,. / DiAt\ ,, D 2 At 

-D 2 At 



4._f,,n+l 4. ,,n+l 4-ii n + l , n+1 A 

~ A V. u i-lj+l,k ^ u ij+l,k ^ u i-lj+2,k ^ u ij+2,kJ 



2h 2 



+ U n+1 flHl^ 



u i-lj-l,k 2 (u i _ 1 j_ lk + u i—1 j +1 k J 



1 1D 2 

= 2 Dl ^ U ^ + 2h| 

+ -f 

4 V +U lj+l,k + u r-lj+2,k + u i!j+2,k 



A i-lj+l,k 



U 



DiAtA n+1 / DjAt 

i - 1 -i- k V^hTJ + U ^ k I 1 + ~W 



n+1 



+ U 



D 2 At 

| - 1 J +1 - k \ 2h 2 
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+u 



-D 2 At 



iJ+2 M 4h 2 



+ u 



-D 2 At 



+u 



n+l 

i-l,]+2,ky 4h 2 



i-l,j+l.k ^ 4h 2 

-D 2 At 



+ u 



-D 2 At 



+ u 



iJ+1 ' k V 4h 2 
-DiAt 



1 1D 2 



i+1 -i'H 2h 2 
3 



-u 



i-l,j,k 



■M-l,i+l,k 



+ a l u ij+l,k + u i-lj+2,k + u ij+2,kJ 



+ U "i,k 



f 



ij,k 



\ 



D,At D,At 
1+— '— + ^— 
K 2h, 



f \ 

D.At 



4h; 



\ 

D.At 



■lj,k 



2/r 

V ~b~ J 



ij + 2,k 



-D.At 



+u 



4h; 
■ ~T~ j 

1 D 



-D t At 

2h\ 

V b J 



1 ^.2 n 
~ D l 5 x U iik + 

2 •'■ 2 h 2 



in n \ ~ n ( n n \ 

V 'J.k i,l-2.l ) ij.k V 'J+2,t 1,7,*/ 



i,7,Jr 



rf "/,* 



<=> au i,i-2,k + bu i-ij,k + cu iJik + bu i+1Jik + au iJ+2ik = d? jjk 
This procedure leads to a linear system Au = f , where 

u = [ u 2,2 U 32 • . . U,V2-1,2 i U23U33 . . . Un/2-1,3 , . . . 

U2,n/2-lU3, n /2-l • • • Un/2-l,n/2-l ]T 

f = [d2,2 d3,2 . . • d,v2-l,2 ! d 2 ,3 (I33 . . . d,v2-l,3 , • ■ ■ d2.11/2-l 
d3,n/2-l ■ • • dn/2-l,n/2-l ]T 

The index £ + 1 is omitted in the vectors u,/for the sake of 
better presentation. The matrix A is given by 

cb a 



b c b 








a 


06 c 


b 






a 


6 


c 


h 









b 


cb 






ii 




be 


b 




ii 




b 


c 


b 


a 






6 


c b 


ll a 








b c b 



b c b 



a b c b 

qOOOOO&c 



For the (i + l)st step of the iterative operator splitting 
method, as restriction operator we choose the most obvious 
restriction operator, the injection. We prefer it instead of the 
standard full weighting operator because it leads to a linear 
system with better structure that is easier to handle. Injection 
is defined by 



„2h _ v 2h 

v ij - v 2i,2j- 



Step i+1: 



,,11+1 _ ,,11 

u i+lj,k+l u i+lj,k+l 



At 



_ 1 
~ 2 

+ D, 



RD, 



,n+l 
i-l,i,k 



n+l , ,,n+l 



2u^+u 



i+l,i,k 



u 



n+l 



i+lj-l,k+l 



2u 



h 2 

"x 

n+l 



+ U 



n + l 



i+lj,k+l T u i+lj+l,k+l 



h 2 



+ ( RDi 



-4-1,1 ,k 



2u u,k + u r+i,i 



h 2 



+ D 2 



_ 9, ,n _i_ „n 

i+l,j-l,k+l z,u i+l,j,k+l "r u i+l,j+l,k+l 

ll 2 



= -^RSjuf+i +-^65uf J t l , VJ . 1 + ^R5 2 u" 



-— x u ij,k 
2 ,J Y u i+lj,k+l 

D 2 At 



D 2 , 
+ -f5 2 u n 



__1 

2 J y u i+ij,k+i "•" 2 



x u i,j,k 



,,n+l I *■ 1 , D,,n + l 

u i+lj-l,k+l I 2h 2 / i-L. 



,n+l ( i + ,,n + l 

j,k \ ov,2 + u i+lj,k+l 



2h 2 



L1 D,At ,, D 2 At 

I r),,n+l 1 | ,,n+l L 

+ Ku ij,k h2 + u i+lj,k+l h 2 
"x u y 



+ u!+ +1 



i,i+i,k+i^ 2h 2 



-D 2 At 



+ Ru 



-DiAt 



i+1 <i' k l 2h 2 



= -(D 1 R5 2 u^ k + D 2 5 2 uf +1Jk+1 ) 



,n+l 



+ U i+lj,k+l 

-D 2 At 



+ u" 



-DjAt 



-l+lj-ljc+l y 2h 2 J - "i-l,],k+l V 2h 2 

± , DiAt 



+ u 



n+l 
i+lj,k+l 



D 2 At 



ij,k+l h 2 ^ u i+lj,k+l h 2 
"x J1 y 



+ U 



n+l 



i+lj+l,k+l 1 2 n 2 



-D 2 At 



+ u 



n+l 



-D^t 



f N 

-D.At 



ij,k+l 



i+uk+iy 2h 2 



= 2( D i 6 x u "i,k+i + D 25>r+i,j, k +i) 

-I- ll n 

~ u i+l,j,k+l 



2 A 



;,7,t+i 



1+ 



DjAf ZJ^Af 



\ 



y+ijk+i 



2h 



v 6 y 



/ -2 n _ ,2 n \ n 

D,& u. . , +D,5 u. , .. , l+u. . , 

V 1 X i,J,k 2 y i+l,j,k+l / l.J.k 



df 



J Ml 



2/i 

X 

au ij-l,k+l + t> u i-lj,k+l + cu ij,k+l + t> u i+lj,k+l + 
au i,j+l,k+l = dij.k+l 

Au = f , where 



u — t u 2,2 U3,2 . • . Un/2-1,2 , U2,3 U3,3 • • ■ Ud/2-1,3 > ■ • • U 2j n/2-l 
U3,n/2-l ■ ■ ■ Un/2-l,n/2-l ]T 

f = [d2,2 d 3?2 . . . d„/2-l,2 , d 2 ,3 d 3j3 . . . dn/2-1,3 , ■ ■ ■ d2,n/2-l 
d3,n/2-l ■ • • dn/2-l,n/2-l ]T 

(The index k+\ is omitted in the vectors u, /for the sake 
of better presentation.) 
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cbQ 

bcb 

06c b 

b c 

a b 

a 



a ii 
n 

6 

c b 
b c 



bcb 



a h r :■' 

a be 



» an n+1 4- hn n+1 4- rn n+1 4- an n+1 4- hn n+1 

^ du ij-l,k + DU i-lj,k-l + cu ij,k + au ij+l,k + DU i+lj,k+l 



,n+l 

-l.i 

= d n 



i.i 



The restriction operators are given by 

Ru» +1 lead to -M^Lk-i + 2u^tU-i + < 



+i 

+l,j-l,k-l 



+ 



9„n + l _i_ a,, n+1 _i_ 9.. n+1 . „n+l . 

^ u i-l,j,k-l "i" ^ u ij,k-l "r ^ u i+lj,k-l "i" u i-lj+l,k-l ~r 



2u 



n + 1 
i,j+l,k-l 



+ U 



n+1 I 

i+lj+l,k-lj 



according to the interpretation of the stencil 



1 
16 



12 1 
2 4 2 
12 1 



The Jacobi iterations for the heat equation are given by 

-u (n+1) + 2u (n+1) - u (n+1) 
u[" +1) = At Dl VlJ + 7 2 ^- 

-u* n+1) + 2u (n+1) - u (n+1) 

h 2 l| 

5. Implementation of the Iterative Splitting Method with 
Embedded Two-grid Method 

The implementation for the heat equation is given by 



,,11+1 _ ,,n i 

u i,j,k u i,j,k _ -J- 

At ' ~ 2 



,n+l 



,n+l _i_ ,,n+l 



u i-l,j,k ^ u ij,k ~ u i+lj,k 



+ PD, 



u 



n+1 



i,j-l,k 



« 



. 1 -2u^ i _ 1 +uP 



■ n+1 _i_ ,,n+l 

-Mj+l,k-l 



+ D t R 



h 2 

"y 

u i-lj,k — 2u L j k + U i+1 j k 



hi 



+ D 



U i,j-l,k+l ^Uji k+1 + Ujj +1 _ k+1 



h 2 



- ^ 6 x u ij,k + ^^ d y u ij,k-l + ^T K6 x u ij,k 
D 2 , 

+ ySXj,k+i 

/-D 2 At\ ,, /-DiAt\ 

~ P <-Wi hhV + Ru ^' k hH + <* 



+ u 



DiAt 



+ Pu 



n+1 2 



D, At 



n+1 _ 

ij,k h 2 n ru ij,k-l h 2 



+ Pu; 



n + 1 



+ U 



ij+l,k-l^ 2h 2 

n+1 

ij,k 



D 2 At\ ,, /-DiAt 

" * + Ru" +1 ' ] 



+1 -i- k V 2h 2 



= o(DiR5^",k + D 2 P5>",k+i) + < 



f \ 

DM 



ij-U 



2/r 



D,At 



i-U*-l 



+11 



V a y 



W+ 1 .* 



-D 2 A* 

2/r 
V a y 



+ L1 



2/i 
V ~7 . 

D,Af 



j,i+l 



/ \ 

D,At ZXAf 

1+— — + — — 



+lj,*-l 



2/r 
v "7" y 



" (A^Uk +Z) 2 8 v U i n + l,j,k + l ) +"",,* : + l 



V. NUMERICAL EXPERIMENTS 

In the following examples, we deal with different test 
examples and their underlying multi-scale physics. 

A. Simple Heat Equation 

We deal with a PDE which is parabolic and has stiffness in 
the diffusion part. 

We have the following equation: 

d t u 1 = D 11 3 xx u 1 + D 21 d xx u 2 , 

3 t u 2 = D 12 a xx u 1 + D 22 a xx u 2 , (80) 

u 1 (0) = u 10 ,u 2 (0) = u 20 (initial conditions), (81) 

% (x, t) = u 2 (x, t) = (boundary conditions), (82) 

where Dll , D21 , D12 , D22 e IR+ and Dll , D12 < 
D21, D22 are the diffusion operators. 

We apply the finite difference scheme for the spatial 
derivatives: 

5xx = [l-2 l]anddyy=[l -2 l]t. 

With the standard pro jection and restriction: 

1 2 11 
R H = 1/16 



2 4 2 
12 1 



P h = 4R H 

The time-derivations are discretised by implicit Euler 
methods and we obtain the following system of linear 
equations: 



Ccf 



Ci n )/At = ACi n+1 + PBcfJi\with q(t n ) = c n (83) 



(cP/i 1 - q+i)/At = RA Cl n+1 + Bcf+Y, 

withCj + l(t n ) = c n 



(84) 



We obtain, after the insertion of the operators, a system of 
linear equations: 



AU 



n+1 _ 



BUf_Y + cup, 



(85) 



This is solved by a two-grid method. Here we have two 
scales and decouple: 



d t u = Au + Bu, 
u(0) = (u 10 ,u 20 ) T , 



(86) 
(87) 
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where u(t) = (uj (t), uj (t)) 
Our split operators are 



T 



for t 6 [0, T]. 



A = fDllflxxO^ R _ fODziflxxA 
VDl 2 3xx0y' V0D 2 2 3xxy 



(88) 



We chose such an example to have AB =£ BA, and 
therefore we have a splitting error of first order for the usual 
sequential splitting methods, called A-B splitting. 

Our numerical results based on two-grid methods are 
presented in the following Table I . 

We choose Dll = D12 = 0.5 and D21 = D22 = 0.05 on the 
time interval [0, 1]. As second order method we choose 
Crank-Nicolson with 6 = 1/2. As fourth order method we 
choose the Gauss Runge-Kutta. We apply a two-grid method, 
see subsection IV-A. 

The numerical results are given in Table I . We apply the 
Z/2 error based on the exact and the numerical solution. 

TABLE I NUMERICAL RESULTS FOR THE FIRST EXAMPLE WITH THE ITERATIVE 
SPLLTTING METHOD AND SECOND AND FOURTH ORDER METHODS 



Number 
of time- 
pajtitions 


Iterative 
Steps 


erri 
(order 2) 


err2 
(order 2) 


erri 
(order 4) 


(order 4) 


2 
2 
2 


1 

10 
100 


4.532le-002 
3.9664e-003 
3;9204e-004 


3,6077e-003 
4.7396e-004 
4.8078e-005 


4.532le-O02 
3-966*5-003 

3.9204e-004 


3.6077e-003 
4.7397e-004 
4.8083e-005 


3 
3 
3 


1 

10 
100 


4.6126e-O04 
7.8129e-006 
8.5988e-008 


3.6077e-003 
2.9285e-005 
2.8270e-007 


4.612Se-004 
7.8O69e-0O6 
8.0050c-008 


3.6077e-003 
2.9289e-005 
2.8682e-007 


4 
4 
4 


1 

10 
100 


4,6126e-004 
4.1883e-007 

5.9521e-O09 


2.2459e-005 
4.2629e-008 
5.48466-009 


4.6126e-004 
4.1321e-007 
4.0839e-010 


2.2464e-005 
4.8154e-008 
4.9968e-011 



Remark 4. In the experiments, we improved the iterative 
splitting scheme with higher order discretisation methods. 
Further we can accelerate the splitting scheme with an 
underlying two-grid method that is embedded in the splitting 
scheme. Both higher order time-discretisation and also 
two-grid methods are necessary to obtain such results. 

B. Second Example: Transport— Reaction Equation 
We deal with the following system of transport equations: 

d t u t + v 1 a x U! = -}U! + \iu 2 , (89) 

d t u 2 + v 2 d x u 2 = \iu x - }u 2 , (90) 

where we have Dirichlet boundary conditions and U\$ , 1*2,0 
are the initial conditions. 

We split the operator into two operators: 

A = ( - Vl _ fl v x 2flx ) - B = (-_£) (91) 

We see that for (i ~ the operators commute. 

We choose vl = 1, v2 = 0.5, 1 = 0.01. We let t e [0, 40] and* 

6 [0,40]. 

We set 

1 4 

At = — and Ax = — 
25 10 

In Fig. 2, we choose u = 0.001 and see that we could say u 
~ and obtain nearly the same results as for the A-B splitting. 



In Fig. 3, we choose u = 0.01 and see that u = and obtain 
more optimal results for the iterative schemes. 

The result is given in Figs. 2-4. 





Fig. 2 Blue - Iterative Operator Splitting (4 iterations); Green - A-B Splitting 




Fig. 3 T=2: Blue - Iterative Operator Splitting (4 iterations); Green - A— B 
Splitting 




Fig. 4 T=60: Blue - Iterative Operator Splitting (4 iterations); Green - A-B 
Splitting 

Remark 5. In this example, we concentrate on the 
comparison between the standard A-B splitting and the 
iterative splitting methods. We obtain more accurate results 
for non-commuting problems. Such problems are related to 
multi-scale problems and we can apply our embedded 
multi-grid method. We have the benefit of receiving higher 
order results without reducing the time-step. 

VI. CONCLUSIONS AND DISCUSSION 
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We presented iterative splitting methods with embedded 
multi-grid and multi-stepping methods. Here the idea was to 
achieve more accurate and faster results than the standard 
splitting schemes (e.g., the A-B splitting schemes). We 
obtained a benefit from this embedded method in accelerating 
the iterative schemes and achieved much more accurate results. 
In the future we will concentrate on nonlinear and matrix 
dependent splitting schemes. 
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